rm(list=ls())
host = "Applejack" #"Ridgeside"
if (host == "AppleJack") {
 setwd("~/github/bmc_netwk_aging_manuscript/R1/1.kaeberlein04plos")
}
if (host == "Ridgeside") {

}

library('flexsurv')
## Loading required package: survival
library('stringr')
source("../lifespan.r")
list.files()
##  [1] "0.example.histogram.density.overlay.Rmd"
##  [2] "0.fit.kaeberlein04.html"                
##  [3] "0.fit.kaeberlein04.pdf"                 
##  [4] "0.fit.kaeberlein04.Rmd"                 
##  [5] "0.histogram.density.overlay.Rmd"        
##  [6] "092304.merged.rls.csv"                  
##  [7] "1.boostrap.kaeblein04.html"             
##  [8] "1.boostrap.kaeblein04.pdf"              
##  [9] "1.boostrap.kaeblein04.Rmd"              
## [10] "2009"                                   
## [11] "bootstrap"                              
## [12] "Kerberlein04PLoS.pdf"                   
## [13] "lifespan.20160926.r"                    
## [14] "lifespan.r"                             
## [15] "sandbox"
tb.ori = read.table("092304.merged.rls.csv", header = T, sep="\t")
#summary(tb.ori)

Explore the fitting outcomes of ‘flexsurv’.

strains = names(tb.ori)
strains
##  [1] "X031104.BY2N.RLS.tab"            "X031104.SK1.RLS.tab"            
##  [3] "X042204.KLBY4716.RLS.tab"        "X052604.sir2D.rls.tab"          
##  [5] "X052604.W303.rls.tab"            "X053104.BY4741.rls.tab"         
##  [7] "X053104.BY4742.rls.tab"          "X053104.JSBY4741.rls.tab"       
##  [9] "X061004.srp1ts.rls.tab"          "fig1a.BY4742"                   
## [11] "fig1a.fob1"                      "fig1a.gpa2"                     
## [13] "fig1a.gpr1"                      "fig1a.hxk2"                     
## [15] "fig1b.BY4742"                    "fig1b.fob1"                     
## [17] "fig1b.hxk2"                      "fig1b.fob1.hxk2"                
## [19] "fig1c.BY4742"                    "fig1c.fob1"                     
## [21] "fig1c.gpa2"                      "fig1c.fob1.gpa2"                
## [23] "fig2a.BY4742"                    "fig2a.sir2"                     
## [25] "fig2a.hxk2.sir2"                 "fig2a.gpa2.sir2"                
## [27] "fig2b.BY4742"                    "fig2b.sir2"                     
## [29] "fig2b.sir2.fob1"                 "fig2c.BY4742"                   
## [31] "fig2c.sir2.fob1"                 "fig2c.sir2.fob1.hxk2"           
## [33] "fig2d.BY4742"                    "fig2d.sir2.fob1"                
## [35] "fig2d.gpa2.sir2.fob1"            "fig3.by4742.2glucose"           
## [37] "fig3.sir2.fob1.2glucose"         "fig3.by4742.05glucose"          
## [39] "fig3.sir2.fob1.05glucose"        "fig3.WT.01glucose"              
## [41] "fig3.sir2.fob1.01glucose"        "fig3.by4742.005glucose"         
## [43] "fig3.sir2.fob1.005glucose"       "fig4a.PSY316"                   
## [45] "fig4a.PSY316.fob1"               "fig4a.PSY316.sir2ox"            
## [47] "fig4b.BY4742.2glucose"           "fig4b.BY4742.005glucose"        
## [49] "fig4b.by4742.SIR2.ox.2glucose"   "fig4b.by4742.SIR2.ox.005glucose"
my.strains = c("fig4b.BY4742.2glucose",  "fig4b.by4742.SIR2.ox.2glucose",  "fig2a.sir2", "fig2b.sir2",              "fig1b.BY4742", "fig1b.fob1","fig1b.hxk2","fig1b.fob1.hxk2"   )
#my.strains = strains
tb = tb.ori[, my.strains]

report = data.frame(my.strains)
report$samplesize = NA; report$R=NA; report$t0=NA; report$n=NA; report$G=NA; #report$longfilename=NA; 

Now, fit all RLS data sets by strains

for( i in 1:length(report[,1])){
  #report$samplesize[i] = length(tb[,i])
  my.data = tb[,i]
  my.data = my.data[! is.na(my.data)]
  report$samplesize[i] = length(my.data)
  
  GompFlex = flexsurvreg(formula = Surv(tb[,i]) ~ 1, dist = 'gompertz')
  WeibFlex = flexsurvreg(formula = Surv(tb[,i]) ~ 1, dist = 'weibull')

  report$avgLS[i] =  mean(tb[,i], na.rm=T)
  report$stdLS[i] =  sd(tb[,i], na.rm = T)
  report$CV[i] = report$stdLS[i] / report$avgLS[i]

  report$GompGFlex[i] = GompFlex$res[1,1]
  report$GompRFlex[i] = GompFlex$res[2,1]
  report$GompLogLikFlex[i] = round(GompFlex$loglik, 1)
  report$GompAICFlex[i] = round(GompFlex$AIC)

  report$WeibShapeFlex[i] = WeibFlex$res[1,1]
  report$WeibRateFlex[i] = WeibFlex$res[2,1]
  report$WeibLogLikFlex[i] = round(WeibFlex$loglik, 1)  
  report$WeibAICFlex[i] = round(WeibFlex$AIC)

  #set initial values
  Rhat = report$GompRFlex[i]; # 'i' was missing. a bug costed HQ a whole afternoon.
  Ghat = report$GompGFlex[i];
  nhat = 6;  
  t0= (nhat-1)/Ghat;
  fitBinom = optim ( c(Rhat, t0, nhat),  llh.binomialMortality.single.run, 
                     lifespan=tb[,i], 
                     #method='SANN') #SANN needs control  
                     method="L-BFGS-B", 
                     lower=c(1E-10, 1, 1), upper=c(1,200,20) );  
  report[i, c("R", "t0", "n")] = fitBinom$par[1:3]
  report$G[i] = (report$n[i] - 1)/report$t0[i]
}

Show the results

#report[ grep("tBY", report$strains), ]
report
##                      my.strains samplesize           R        t0        n
## 1         fig4b.BY4742.2glucose         60 0.002388227  35.59111 7.703055
## 2 fig4b.by4742.SIR2.ox.2glucose         60 0.003068054  66.82089 8.067874
## 3                    fig2a.sir2         90 0.002656909  16.80162 8.109348
## 4                    fig2b.sir2         90 0.002656909  16.80162 8.109348
## 5                  fig1b.BY4742        250 0.004796452  58.10834 8.109395
## 6                    fig1b.fob1        140 0.003056547  71.37470 7.695271
## 7                    fig1b.hxk2        120 0.005736777 103.59300 7.578456
## 8               fig1b.fob1.hxk2        160 0.005902597 120.65937 6.023612
##            G    avgLS     stdLS        CV  GompGFlex   GompRFlex
## 1 0.18833511 26.06667  7.557389 0.2899254 0.13946781 0.002283890
## 2 0.10577341 34.60000 10.834972 0.3131495 0.07468834 0.004042203
## 3 0.42313480 13.96667  3.491402 0.2499810 0.27970246 0.003555048
## 4 0.42313480 13.96667  3.491402 0.2499810 0.27970246 0.003555048
## 5 0.12234725 26.62400  9.418144 0.3537464 0.08544660 0.006577615
## 6 0.09380453 37.75000 13.405746 0.3551191 0.06987148 0.003424637
## 7 0.06350290 36.73333 16.338197 0.4447785 0.04823602 0.006725196
## 8 0.04163466 48.28125 21.337707 0.4419460 0.04143860 0.004249514
##   GompLogLikFlex GompAICFlex WeibShapeFlex WeibRateFlex WeibLogLikFlex
## 1         -206.8         418      3.925269     28.75636         -206.2
## 2         -235.6         475      3.369961     38.51039         -227.8
## 3         -246.0         496      4.431406     15.29449         -240.4
## 4         -246.0         496      4.431406     15.29449         -240.4
## 5         -937.3        1879      3.024640     29.77138         -914.0
## 6         -565.5        1135      3.139805     42.25834         -559.6
## 7         -507.6        1019      2.397501     41.42400         -502.2
## 8         -713.6        1431      2.441562     54.41078         -714.3
##   WeibAICFlex
## 1         416
## 2         460
## 3         485
## 4         485
## 5        1832
## 6        1123
## 7        1008
## 8        1433

Output

# write.csv(report, file = 'sandbox/_report_kaeberlein04_fit.csv', row.names = FALSE)

Overlay the binomial aging model with observation.

see http://hongqinlab.blogspot.com/2013/12/binomial-mortailty-model.html m = R ( 1 + t/t0)^(n-1) s = exp( (R t0/n)*(1 - (1+t/t0)^n ) )

for( i in 1:length(report[,1])){
  #report$samplesize[i] = length(tb[,i])
  my.data = tb[,i]
  my.data = my.data[! is.na(my.data)]

  ret = calculate.s(my.data)
  plot( ret$s ~ ret$t, main=my.strains[i]);
  print (report[i,  ]);

    #overlay binomial aging viability
  print (report[i, c("R", "t0", "n", "G")] );
  t = seq(0,max(ret$t*1.1), by=0.1);
  # s = exp( (R t0/n)*(1 - (1+t/t0)^n ) )  
  s = exp( (report$R[i]* report$t0[i]/report$n[i])*(1 - (1+t/report$t0[i])^report$n[i] ) ) ; 
  lines(s ~ t, col='red')
  
  #overlay gompertz  viability
  s.g =  G.s( c(report$GompRFlex[i], report$GompGFlex[i], 0), t )
  lines(s.g ~ t, col='blue')

}  

##              my.strains samplesize           R       t0        n         G
## 1 fig4b.BY4742.2glucose         60 0.002388227 35.59111 7.703055 0.1883351
##      avgLS    stdLS        CV GompGFlex  GompRFlex GompLogLikFlex
## 1 26.06667 7.557389 0.2899254 0.1394678 0.00228389         -206.8
##   GompAICFlex WeibShapeFlex WeibRateFlex WeibLogLikFlex WeibAICFlex
## 1         418      3.925269     28.75636         -206.2         416
##             R       t0        n         G
## 1 0.002388227 35.59111 7.703055 0.1883351

##                      my.strains samplesize           R       t0        n
## 2 fig4b.by4742.SIR2.ox.2glucose         60 0.003068054 66.82089 8.067874
##           G avgLS    stdLS        CV  GompGFlex   GompRFlex GompLogLikFlex
## 2 0.1057734  34.6 10.83497 0.3131495 0.07468834 0.004042203         -235.6
##   GompAICFlex WeibShapeFlex WeibRateFlex WeibLogLikFlex WeibAICFlex
## 2         475      3.369961     38.51039         -227.8         460
##             R       t0        n         G
## 2 0.003068054 66.82089 8.067874 0.1057734

##   my.strains samplesize           R       t0        n         G    avgLS
## 3 fig2a.sir2         90 0.002656909 16.80162 8.109348 0.4231348 13.96667
##      stdLS       CV GompGFlex   GompRFlex GompLogLikFlex GompAICFlex
## 3 3.491402 0.249981 0.2797025 0.003555048           -246         496
##   WeibShapeFlex WeibRateFlex WeibLogLikFlex WeibAICFlex
## 3      4.431406     15.29449         -240.4         485
##             R       t0        n         G
## 3 0.002656909 16.80162 8.109348 0.4231348

##   my.strains samplesize           R       t0        n         G    avgLS
## 4 fig2b.sir2         90 0.002656909 16.80162 8.109348 0.4231348 13.96667
##      stdLS       CV GompGFlex   GompRFlex GompLogLikFlex GompAICFlex
## 4 3.491402 0.249981 0.2797025 0.003555048           -246         496
##   WeibShapeFlex WeibRateFlex WeibLogLikFlex WeibAICFlex
## 4      4.431406     15.29449         -240.4         485
##             R       t0        n         G
## 4 0.002656909 16.80162 8.109348 0.4231348

##     my.strains samplesize           R       t0        n         G  avgLS
## 5 fig1b.BY4742        250 0.004796452 58.10834 8.109395 0.1223473 26.624
##      stdLS        CV GompGFlex   GompRFlex GompLogLikFlex GompAICFlex
## 5 9.418144 0.3537464 0.0854466 0.006577615         -937.3        1879
##   WeibShapeFlex WeibRateFlex WeibLogLikFlex WeibAICFlex
## 5       3.02464     29.77138           -914        1832
##             R       t0        n         G
## 5 0.004796452 58.10834 8.109395 0.1223473

##   my.strains samplesize           R      t0        n          G avgLS
## 6 fig1b.fob1        140 0.003056547 71.3747 7.695271 0.09380453 37.75
##      stdLS        CV  GompGFlex   GompRFlex GompLogLikFlex GompAICFlex
## 6 13.40575 0.3551191 0.06987148 0.003424637         -565.5        1135
##   WeibShapeFlex WeibRateFlex WeibLogLikFlex WeibAICFlex
## 6      3.139805     42.25834         -559.6        1123
##             R      t0        n          G
## 6 0.003056547 71.3747 7.695271 0.09380453

##   my.strains samplesize           R      t0        n         G    avgLS
## 7 fig1b.hxk2        120 0.005736777 103.593 7.578456 0.0635029 36.73333
##     stdLS        CV  GompGFlex   GompRFlex GompLogLikFlex GompAICFlex
## 7 16.3382 0.4447785 0.04823602 0.006725196         -507.6        1019
##   WeibShapeFlex WeibRateFlex WeibLogLikFlex WeibAICFlex
## 7      2.397501       41.424         -502.2        1008
##             R      t0        n         G
## 7 0.005736777 103.593 7.578456 0.0635029

##        my.strains samplesize           R       t0        n          G
## 8 fig1b.fob1.hxk2        160 0.005902597 120.6594 6.023612 0.04163466
##      avgLS    stdLS       CV GompGFlex   GompRFlex GompLogLikFlex
## 8 48.28125 21.33771 0.441946 0.0414386 0.004249514         -713.6
##   GompAICFlex WeibShapeFlex WeibRateFlex WeibLogLikFlex WeibAICFlex
## 8        1431      2.441562     54.41078         -714.3        1433
##             R       t0        n          G
## 8 0.005902597 120.6594 6.023612 0.04163466

overlay probility mass function with binomial model

see http://hongqinlab.blogspot.com/2013/12/binomial-mortailty-model.html m = R ( 1 + t/t0)^(n-1) s = exp( (R t0/n)(1 - (1+t/t0)^n ) )
pdf = s
m

for( i in 1:length(report[,1])){
  #report$samplesize[i] = length(tb[,i])
  my.data = tb[,i]
  my.data = my.data[! is.na(my.data)]
  
  h= hist(my.data, br=max(my.data)/2)
  plot( h$density ~ h$mids, main=my.strains[i], xlab="RLS",ylab="density")
  t= seq(0, max(h$mids), by=0.1)
  s = exp( (report$R[i]* report$t0[i]/report$n[i])*(1 - (1+t/report$t0[i])^report$n[i] ) ) ; 
  m = report$R[i]*(1 + t/ report$t0[i])^(report$n[i] -1 )
  pdf = s*m
  lines( pdf ~ t, col='red')
  
  s.g =  G.s( c(report$GompRFlex[i], report$GompGFlex[i]), t ); 
  m.g =  report$GompRFlex[i]*exp(report$GompGFlex[i]*t)
  pdf.g = s.g * m.g
  lines( pdf.g ~ t, col="blue")
}

It seems binomial model aging is a reasonable fit whenevern Gompertz model is a reasonable fit.

for( i in 1:length(report[,1])){
  my.data = tb[,i]
  my.data = my.data[! is.na(my.data)]
  
  h= hist(my.data, br=max(my.data)/2);
  hist(my.data, probability = TRUE, col='black', border='white', xlab='RLS', ylab='probability',
       main=my.strains[i], ylim=c(0, max(h$density)*1.1), xlim=c(0, max(h$mids)*1.1))
  #plot( h$density ~ h$mids, main=my.strains[i], xlab="RLS",ylab="density")
  #par(new=TRUE);
  t= seq(0, max(h$mids)*1.2, by=0.1)
  s = exp( (report$R[i]* report$t0[i]/report$n[i])*(1 - (1+t/report$t0[i])^report$n[i] ) ) ; 
  m = report$R[i]*(1 + t/ report$t0[i])^(report$n[i] -1 )
  pdf = s*m
  lines( pdf ~ t, col='red')
  
  s.g =  G.s( c(report$GompRFlex[i], report$GompGFlex[i]), t ); 
  m.g =  report$GompRFlex[i]*exp(report$GompGFlex[i]*t)
  pdf.g = s.g * m.g
  lines( pdf.g ~ t, col="blue")
}